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ABSTRACT 

We carried out 2D-axisymmetric MHD simulations of core-collapse supernovae for rapidly-rotating 
magnetized progenitors. By changing both the strength of the magnetic field and the spatial resolution, 
the evolution of the magnetorotational instability (MRI) and its impacts upon the dynamics are 
investigated. We found that the MRI greatly amplifies the seed magnetic fields in the regime where 
not the Alfven mode but the buoyant mode plays a primary role in the exponential growth phase. 
The MRI indeed has a powerful impact on the supernova dynamics. It makes the shock expansion 
faster and the explosion more energetic, with some models being accompanied by the collimated-jet 
formations. These effects, however, are not made by the magnetic pressure except for the collimated- 
jet formations. The angular momentum transfer induced by the MRI causes the expansion of the 
heating region, by which the accreting matter gain an additional time to be heated by neutrinos. The 
MRI also drifts low-1^ matter from the deep inside of the core to the heating region, which makes the 
net neutrino heating rate larger by the reduction of the cooling due to the electron capture. These 
two effects enhance the efficiency of the neutrino heating, which is found to be the key to boost the 
explosion. Indeed we found that our models explode far more weakly when the net neutrino heating is 
switched off. The contribution of the neutrino heating to the explosion energy could reach 60% even 
in the case of strongest magnetic field in the current simulations. 

Subject headings: supernovae: general — magnetohydrodynamics (MHD) — Instabilities — methods: 
numerical — stars: magnetars 


1. INTRODUCTION 

Magnetic field and rotation are ubiquitous in stars. 
MiMeS survey has observed over 550 Galactic O- and 
B-type stars, and detected the surface magnetic fields of 
> 100 G for ^ 10 % of them (see Wade & the MiMeS Col¬ 


laboration (2014) for review). Estimating the upper limit 
of the currently-undetected magnetic fiel d to be ^ 100 G, 
Wade & the MiMeS Collaboration! (|2014|) argued that the 
distribution of the magnetic fields tor uiassive stars may 
be bimodal: a small population of strong magnetic fields 
(> 1 kG) and a large majority of weak magnetic fields 
(< 100 G). A magnetic field of 1 kG corresponds to the 
magnetic flux of ^ 10^^ G cm^ f or a 17 Mq star with the 
radius of ^ 8Rq (McNally|1965 ), which is comparable to 
t hat of magnetars. 


Ramfrez-Agudelo et al. (2013) measured the surface 
rotational velocities of 216 0-type stars, and found that 
25 % of the sample have vsini > 200 km s“^ while the 
rest of them are slow rotators. According to st ellar evo¬ 
lution calculations by Woosley & Heger (2006), if a star 
is rotating fast enough initially, the rotational mixing 
prevents a very efficient angular momentum transport 
between the helium core and the hydrogen envelope, and 
the central iron core maintains a large amount of an¬ 
gular momentum at pre-collapse stage. They inferred 
the rotation period of a neutron star to be 2.3-9.7 ms 
for such evolutions of a 16 Mq star with solar metallic- 
ity. Then, the high-rotational-vel ocity population found 
by Ramfrez-Agudelo et al. (2013) might produce proto- 
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neutron stars rotating with a period similar to those of 
millisecond pulsars (MSPs). 

The influences of magnetic field and rotation on core¬ 
collapse supernovae have been studied as a possible agent 
to drive explosion other than neutrino heating, while the 
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). MHD 

core-collapse simulatic 
main focus on rather ( 
10^^ G and ^ ^ 

correspond to the m^ 
MSP-class rotation (e. 

)ns done so far have placed the 
extreme cases, viz., ^ 10^^- 

L rad s“^ at pre-collapse, which 
rgnetar-class magnetic field and 
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mat ions, tne mag- 
netic field wound by differential rotation grows to dy¬ 
namically important strengths and later drives a strong 
outflow along the rotation axis, reproducing the typical 


supernova-explosion energy of 


exp 


10^^ erg. 


Since the magnetic field and rotation in massive stars 
are likely to have wide range of values as mentioned 
above, it may be also important to study more ’’ordi¬ 
nary” cases. In these cases amplification mechanisms 
that are more efficient than the simple winding are im¬ 
perative to produce the field strength of 10^^ G out¬ 
side the proto-neutron star, which may be necessary to 
imp act on the supernova dyna mics. For non-rotating 


case 


Endeve et al.| ( |2Q1Q[ |2Q12[ ) numerica lly studied the 


stand i ng acc retion shock instability, while [Obergaulinger| 


et al. (2014) investigated the convection. In both cases, 
the aniplification is rather modest, and the impacts on 
dynamics are found to be minor. 

If the iron core is initially rotating rapidly, another 
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candidate of an efficient field amplification mechanism 
in core-collapse supernovae is the magnetorotational in¬ 
stability (MRI), which basically occurs in differentially 
rotat ing systems ( [Balbns fc Hawley|1991 Akiyama et ah 
2003). Simulations of the MKi for weak seed magnetic 
fields are computationally demanding, since the wave¬ 
length of the fastest growing mode is quite small com¬ 
pared with the size of the iron core, ^ 1000 km: 
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where va is Alfven velocit30 In fact, most of pre¬ 
vious core-collapse simulations assuming sub-magnetar- 
class magnetic fields have insufficient spatia l resolutions 


to captu re the MRI dMoiseenko et ah 2006 


Burrows et 


^[2007 Takiwaki et al.||2009| ]^ In order to resolve the 


fastest growing mode, local si mulation boxes are utilized 
in some 2D/3p cornputations ([Oberganlinger et al.j2009 


global dynamical, Sawai & Yamada (2014) carried out 


similar but longerAerm simulations up to several hun¬ 
dred milliseconds after bounce, employing the simple 
light bulb approximation for neutrino transfer. They 
found that the MRI indirectly enhances the neutrino 
heating, and thus boost the explosion. Per forming 3D 
sim ulatio ns for a thin layer on the equator, [Masada H 
al. (2015) argued another possible effect of the MRI, i.e., 
the enhancement of neutrino luminosity by MRI-driven 
turbulence around the pr oto-neutron star sur face. 


This pape r is a se quel to Sawai et al. (2013b) and Sawai 


In spite that Moiseenko et al.j (|2L)Uti|> found an exponential 


sion and conclusion are given in Section]^ 


2. NUMERICAL METHODS 


We adopt a ISMq star (jWoosley & Weaver jjl 995 j) for 
the progenitor of core-collapse simulations, adding mag- 
netic fields and rotations by hand. The following ideal 
MHD equations and the equation of electron number 
density are numerically solve d by a time-explicit Eule- 
rian MHD code, Yamazakura ( Sawai et al.||2013a ): 
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Masada et al.||20f^ Guilet Mnller||20f^ Rembiasz et al 

2015p. The problems in the local simulations, however, 
are the difficulties in taking into account the effects from 
a nd feedbacks t o dyna mically changing structures. 


Sawai et al.j (2013b) conducted the first global core- 
collapse simulations for sub-magnetar-class magnetic 
fields with a sufficient spacial resolution to capture the 
MRI albeit in 2D axisymmetry, and found that the mag¬ 
netic field is amplified by the MRI to dynamically im¬ 
portant strengths. In order to study ffs ini pacts on the 


& Yamada (2014). We conducted 2D-axisymmetric 
high-resolution simulations of core-collapse for rapidly- 
rotating magnetized progenitors, changing the initial 
magnetic field strength and the spatial resolution. The 
initial magnetic field strength assumed here, ^ 

10^^ G, are one or two orders of magnitude smaller than 
the extreme values adopted in some previous simulations 
mentioned above. 

The rest of the paper is organized as follows. In Sec¬ 
tion we describe the numerical method and models. 
The results are presented in Section and the discus- 

^ The wavelength of the fastest growing mode given here is one 
obtained for cylindrical rotation laws, with neglecting buoy¬ 

ancy. We deal with the general rotatio n la ws, D(ti7,2;), taking the 
buoyancy into acc ount later in Secti o n |3.I| 


growth of magnetic held and claimed tliat the growth is due to 
a Tayler-type ’’magnetorota tional instability”, which i s completely 
different from one found byjBalbus & Hawleyj (|I99I|). Note how¬ 
ever, that the property of tJie instability is still unclear, and no 
other groups succeeded to reproduce their results to date. 
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where and are the changes of energy den- 

sity due to neutrino/anti-neutrino absorptions and emis¬ 
sions, respectively, and and are the simi¬ 

lar notations for the changes of electron number den¬ 
sity. The other symbols have their usual meanings. 
The electron fraction, Yp, i s give n by the prescription 
suggested by [Liebendorferj (|2005|) until bounce. After 
that, where Liebendorfer's prescription is no longer valid. 
Equation is solved to obtain Yg = where 

rriu = 1.66 X 10“^^ g is the atomic mass unit. We as¬ 
sume Newtonian mono-pole g ravity. A tabulated n uclear 
equation of state produced by Shen et al. ( 1998a|b) is uti¬ 
lized. Gomputations are done with polar coorcfinates in 
two dimensions, assuming axisymmetry and equatorial 
symmetry. 

We take into account interactions of electron neutri¬ 
nos Ee and anti-neutrinos Pg with nucleons. Instead of 
dealing with detailed neutr ino transport, t he lig h t bulb 
approximation is used as in Murphy e t al. (2009); Nord- 
haus et al. (2010); Hanke et al. (|2012|). Taking the ultra- 
relativistic limit for electrons and positrons, assuming 
the Eermi-Dirac distribution with vanishing chemical po¬ 
tential for Ee and Pg, and neglecting the phase space 
blocking, we evaluate the source term related to z/g/pg 
absorption (z^g -h n ^ e“ + p, z/g -h p ^ e+ -h n) in the 
energy equation @ as 
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(Janka 2001), where a = 1.26 is the charged-current 
axial-vector coupling constant, ao = 1.76 x 10 cm^ the 
characteristic cross section of weak interaction, (e^^) = 

20.S the mean square neutrino energy, mgC^ = 
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0.511 MeV the rest-mass energy of electron, the neu¬ 
trino luminosity, r the distance from the center, and 
{fly) the so-called flux factor. Here, we assume the 
same luminosity and spectral temperature for z^e and Pe- 
In the present simulations, Ly^ = 1.0 x 10^^ erg s“^, 
kTy^ =4.0 MeV, and {jiy) = 1.0 are chosen. Similarly, 
the source term related to z^e/^e absorption in the equa¬ 
tion of Tie is 


/<^abs _ 

Vat — 


3a‘^ + 1 (To{eyJ p 


{nieC^Y '^u 47rr^(/ijy) 


(Vn-V^), (8) 


where {ey^) =4.11 (kTy^) is the mean neutrino energy. 

The source term related to the z^e/^e emission (e“ + 
p ^ + e+ + n ^ Pe +p) in Equation Q is given by 
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(|Janka|2001j) , where pe is the electron chemical potential 
normalized by the temperature, and 
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Similarly, the emission source term in Equation ^ 
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The Fermi integrals (Equation ([l0|)) are calculated as 
follows: 
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The source terms given by Equations Q , (§, 0,0 
are valid only in optically thin regions, and must decrease 
toward the optically thick regions. To mimic such reduc¬ 


tio n they are multiplied by e following Murphy et 


al. (2009). Here, the effective optical depth is dehned as 

poo 
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where the effective opacity is given as 
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lO^^g cm“^ J V 4MeV 
from the Equations (10), (11), and (14) of Janka ( | 2001 ) 


Yp), 

(18) 


Before conducting high-resolution simulations to cap¬ 
ture the MRI, we first follow the collapse of the 15 Mq 
progenitors until several 100 ms after bounce by low- 
resolution simulations, whose numerical domain spans 
from the radius of 100 m to 4000 km. We refer to these 
simulations as background (BG) runs. In the BG runs, 
the core is covered with x Nq = 720 x 60 numerical 
grids, where the spatial resolution is 0.4-23 km. 

The pre-collapse cores are assumed to be rapidly ro¬ 
tating with the initial angular velocity profile of 

iQ-ri 

where r is the distance from the center of the core. 
The parameters are chosen as vq = 1000 km and 
flo = 2.73 rad s“^, corresponding to a millisecond proto¬ 
neutron star after collapse. The initial rotational energy 
divided by the gravitational binding energy is 2.5 x 10“^. 

We assume that the pre-collapse magnetic fields have 
dipole-like configurations produced by electric currents 
of a 2 D-Gaussian-hke distribution centered at {vu^z) = 
(tI7o,0), 


j^{w,z)=joe 




VJqVJ 

w‘i + tiY 


where (lu, z) are cylindrical coordinates, 
[w — vjqY + ^ = arccos( 2 :/f), and 
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(ISawai et aLir2013a|). Ghanging jo, we perform three BG 


runs with different strengths of magnetic fields, where the 
maximum strengths at pre-collapse, Hp^e, are 5.0x10^^, 
1 . 0 x 10 ^^, and 2 x 10 ^^ G. Hereafter we refer to these BG 
runs as B5el0bg, Blellbg, and B2ellbg, respectively. 
The rest of parameters are set as zuq = rdec = 1000 km 
and e = 0.5 in all the computations. The initial mag¬ 
netic energy divided by the gravitational binding energy 
is quite small, 2.1 x 10 “^, even for the strongest-held 
model, B2ellbg. We also computed models without 
magnetic held and rotation as well as a model having 
rotation alone for comparison. 

In order to capture the growth of MRI we conduct 
high-resolution simulations with the numerical domain 
spanning 50 < (r/km) < 500 (referred to as MRI runs). 
The initial conditions of MRI runs are given by map¬ 
ping the data of the BG runs onto the above domain at 
5 ms after bounce. In order to satisfy the divergence free 
constraint on the magnetic held, not the mag netic held 
itse lf but t he vector potential is mapped as in [Sawai 


The inner and outer radial boundary condl 
the MRI runs are given by the data of the basic 


al. (|2013al) 
tions lor ti 

runs, except that the inner boundary conditions of are 
determined to satisfy the divergence-free condition. The 
grid spacing is such that the radial and angular grid sizes 
are the same, viz. Ar = rAO, at the innermost and out¬ 
ermost cells. Eor each BG run, four MRI runs with dif¬ 
ferent grid resolutions are carried out. Our choice of the 
resolution at r = 50 km, A 50 , (and the numbers of grids. 
Nr X Ne), is 12.5 m (9250 x 6400), 25 m (4650 x 3200), 
50 m (2300 x 1600), and 100 m (1160 x 800). We label 
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the MRI runs by the initial field strength of the corre¬ 
sponding BG run followed by the spatial resolution. For 
example, the MRI run using the data of model B5el0bg 
and A 50 = 12.5 m is referred to as model B5elOA12.5. 
For a set of models involving the same initial magnetic 
field, we use a term “model series”, e.g., models series 
B5el0. 

In dealing with the MHD equations in the polar coordi¬ 
nates, we should be cautious about numerical treatments 
of the coordinate singularities at the center of the core 
(r = 0) and the pole (^ = 0). In the vicinity of the pole, 
the regularity conditions demand that the expansions of 
ve^ and with respect to 0 should not contain 

6 >-independent terms, which is not necessarily satisfied 
in numerical simulations. In order to numerically meet 
the regularity conditions in the vicinity of the pole albeit 
approximately, we remove the region of 0 < 0.3° from 
the numerical domain and impose boundary conditions 
based on the regularity conditions except for Bq^ which is 
determined by the divergence-free constraint. To diffuse 
undesirable fluctuations that tend to violate the regular¬ 
ity, we further introduce an artificial resistivity only at 
the cells closest to the pole in the form of 


= - j -, 


( 22 ) 


where a is a dimensionless factor, Umax the local maxi¬ 
mum characteristic speed, A the grid width, and Ib the 
scale height of the magnetic field. The factor a is auto¬ 
matically controlled between 0 . 1 - 10 ^ during the simula¬ 
tions depending on how well the regularity condition is 
satisfied. 

In order to maintain the regularity conditions approxi¬ 
mately around the center in the BG runs, we remove the 
central part within the radius of 100 m from the numer¬ 
ical domain and take a similar remedy. 


3. RESULTS 

3.1. The Growth of MRI 

The stability condition of the axisymmetric MRI for 
general rotation laws is given by 
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Figure 1. Color maps of the dominant modes and growth rate 
for model B5elOA12.5 at tpb = 7 (upper panel) and tpb = 12 
(lower panel). The red and blue colors, respectively, represent the 
locations where buoyant mode and Alfven mode are dominant. 
The growth rate is multiplied by -1 for buoyant-mode-dominant 
regions. The boxes in the upper panel correspond to the plot areas 
of Figure 

mode with the wavenumber of 

^FGM • Va = COsOkl^ ^ ^ -, (27) 


and the growth rate of 


C = {GzBz tan^ Ok - tan Ok + + Ilvo) 

>0, (23) 

where Ok is the angle between the perturbation wavenum¬ 
ber /c and the 2 ;-axis. 
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(26) 


(Balbus 


1995 


Obergaulinger 


and T = d \nP/d\np\sXe 
et al.|[ 2 Qi^ . 

The MRI involves two distinct modes, namely, Alfven 
mode and buoyant mode, where the former appe ars for 
C < 0 , and the latter only emerges for C -h 4 < 0 (Urpin 
1996). Which mode dominates over the other for a hxed 


Ok depends on the value of C (| Obergaulinger et al.|2009|) 
For — 8 < C < 0, the fastest growing mode is the Alfven 


<^FGM = COs 6 >/cO—-—. (28) 

For C < — 8 , the fastest growth occurs with 

^FGM = COS Okl^V C + 4, (29) 

for k^GM ' va = 0^ i.e., it is the buoyant mode. 

Since C depends on Ok^ the dominant mode differs for 
different directions. In order to find the dominant mode 
for a fixed spatial point, we vary Ok numerically in the 
range of [— 7 r /2 : 7 r/ 2 ]. The result is shown in Figure for 
model B5elOA12.5 at the postbounce time of tpb = Tand 
12 ms, where the red and blue colors represent buoyant¬ 
mode- and Alfven-mode-dominant regions, respectively 
and the shades of the colors indicate the growth rat^ 
It is evident that the dominant mode is different from 
location to location, and the regions dominated by the 

^ The resolutions of color maps in this paper are not the same 
as those of simulations, where the former are reduced to decrease 
the size of figures. 
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Figure 2. Cumulative volume fractions having A^mri smaller than 
a given value for the all models. 

buoyant-mode have on average larger growth rates than 
those dominated by the Alfven mode. 

The growth of the Alfven-mode, although slower than 
the buoyant mode, still may have an important effect 
on the magnetic field amplification in the locations of 
its dominance. It is hence important to know how 
well we numerically resolve the fastest-growing Alfven 
mode (TGAM), whose wave number is given by Equa¬ 
tion ([27| ). Figure]^ shows for all the models the cumula¬ 
tive fraction of the volume that has A^mri smaller than 
a given value. Here A^mri is defined at each point to be 
the ratio of the wavelength of FGAM to the grid size 
and is measure of how well FGAM is resolved numeri¬ 
cally. We introduce to characterize the models, a factor 
^ = (-Bi„,max/10“G) (A 5 o/ 100 m)“^ since only the ini- 
tial strength of the magnetic field and spacial resolution 
are different among the current set of the models. In 
fact, similar distributions is obtained for the models hav¬ 
ing the same ^ at early epochs when the ma gnetic field 
is a lmost passive (see Figure [^. According to Shibata et 


al. (2006), A^mri ^ 10 is required to capture the linear 
growth of the Alfven mode. In our weakest-field model 
series B5el0, the volume fractions with A^mri < 10 are 
0.6, 0.27, 0.12, and 0.018 for models AlOO = 0.5), 
A50 = 1), A25 = 2), and A12.5 = 4), re¬ 

spectively. We hence believe that our highest-resolution 
models should be able to capture the linear growth of the 
Alfven mode well. 

In order to confirm that both the buoyant and Alfven 
modes are growing indeed in the regions predicted in Fig¬ 
ure!^ we examine the wavelengths of the growing modes 
in these regions. The upper panels of Figure show 
the color maps of the poloidal magnetic field strength 
a-t tpb = 9.5 ms in a region around the equator, indi¬ 
cated by the large box in the upper panel of Figure 
where two red belts of buoyant-mode dominance are ob- 
servecQ We found that the patterns of strong-magnetic- 
field filaments seen in the upper panels of Figure have 
grown in the regions of buoyant-mode dominance: in the 
case of model B5elOA12.5, this region corresponds to 
the right-side red belt observed in the large box indi¬ 
cated in the upper panel of Figure which has been 
advected leftward during tpb = 7.0-9.5 ms. We compare 


® Although the upper panel of Figure [T] is depicted for model 
B5elOA12.5, its feature is very similar among all the models at 
this point of time. 


models B5elOA12.5 and B2ellA12.5, which have differ¬ 
ent initial field strengths. As expected for the buoyant 
mode, the wavelengths of the growing modes, which are 
evaluated from the sizes of the patterns, are nearly iden¬ 
tical between the two models. The wavelengths of the 
growing modes observed here may reflect the scale of the 
dominant perturbation, which may come from numerical 
noises. To see if this is true, we performed test simu¬ 
lations for the two models in which a perturbation of 
u' = Rq sin(27r2:/Aprt) is given at the beginning of the 
MRI runs, where the amplitude and the wavelength of 
the perturbation is set as 1% and 500 m, respectively. As 
a result, we found that the wavelengths of growing modes 
are shorter than those of the models without perturba¬ 
tion (see panels (c) and (d) of Figure]^, which indicates 
that the observed modes depend on the dominant scale 
of perturbations. 

The lower panels of Figure zoom in the area around 
0 = 35° in the vicinity of the inner boundary, where 
a pocket of buoyant-mode-dominant regions are sur¬ 
rounded by an Alfven-mode-dominant r^ion (see the 
small box in the upper panel of Figure CT. The wave¬ 
length of the growing mode there is shorteiTor the weaker 
initial field. In fact, the widths of protruding mag¬ 
netic flux loops in the panel (e) of Figure 1^ for model 
B5elOA12.5 are about three times smaller than those in 
the panel (f) for model B2ellA12.5. The ratio is close 
to four, the value expected for the Alfven mode. With 
these facts, we believe that our simulations capture both 
the buoyant and Alfven modes correctly. 

Figure plots the time evolutions of the magnetic en¬ 
ergies, integrated over the whole numerical domain for 
all models. Those for the MRI runs are plotted until 
the shock surface reaches the outer boundary. Only in 
model B2ellA12.5, the plots are continued for another 
12 ms after the shock surface passes the outer boundary. 
Since the magnetic energies flowing out of the bound¬ 
ary during this 12 ms are found to be negligible, i.e., 
0.81% and 0.24% of the total poloidal and toroidal mag¬ 
netic energies at the end of the plot, respectively, we do 
not take them into account in the following discussion. 
The exponential growth of the energy of the poloidal 
component, is apparent during the first ^10 ms 

for all the MRI runs. In each model series, the growth 
timescale becomes shorter (or the growth rate is larger) 
for higher resolutions until it converges to ^3-3.5 ms. 
These timescales well match the theoretical prediction 
for the buoyant-mode of ^2000 rad s“^, which is shown 
in the upper panel of Figure This implies that the 
exponential growth is dominated by the buoyant mode. 
Indeed, the comparison between the lower panel of Fig¬ 
ure!^ and the upper panel of Figure indicates the co¬ 
incidence of the locations, where the poloidal magnetic 
field is preferentially amplified, with those of buoyant¬ 
mode dominance at tpb = 12 ms, around the end of the 
exponential growth. From the numerical convergence we 
observed, it is suggested that the high spatial resolution 
is required even for the buoyant mode, in which all wave¬ 
lengths grow at an equal rate. 

After the exponential growth phase ceases, Eb^ con¬ 
tinues to increase gradually until it reaches saturation 
roughly around tpb =210, 270, and 160 ms for model 
series B5el0, Blell, and B2ell, respectively (see Fig¬ 
ure . During this phase the region of strong magnetic 
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Figure 3. Color maps for the strength of the poloidal magnetic field for models B5elOA12.5 (left panels) and B2ellA12.5 (right panels). 
The upper four panels zoom in a part of equatorial region (presented by the large box in the upper panel of Figure while the lower 
two panels that of a middle-latitude region (the small box in the upper panel of Figure [^. Panels (c) and (d) are for models with initial 
perturbations. 


field, say, B > 10^^ G, spreads over a considerable vol¬ 
ume inside the radius of ^ 100 km (see the lower panel 
of Figure]^ for B5elOA12.5). 

As can oe seen from each panel of Figure the sat¬ 
urated values of Eb^ do not converge, which may be 
because the turbulen ce is not yet fully c aptured due to 
numerical diffusivity ( Sawai et aL||2013b ). Nevertheless, 
since model series Blell and B2ell show a trend of con¬ 
vergence, we may be able to estimate the converged val¬ 
ues by fitting the time-averaged Eb^ for the different 
resolution models with suitable functions. Taking the 


time averages over tpb = 270-330 ms and 165-185 ms for 
model series Blell and B2ell, respectively, we fitted the 
results with functions in the form of 

(£^Bp,sat) = 01-02 exp (-Oa/Aso), (30) 

where ai, a 2 , and as are the parameters to be deter¬ 
mined. The values obtained for ai are 1.4 x 10^^ and 
3.1 X 10^^ erg for model series Blell and B2ell, re¬ 
spectively (see Figure]^. It is also found that the sat¬ 
urated values of Eb^ tor the highest-resolution models 
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Figure 4. Evolutions of the magnetic energies of the poloidal and 
toroidal components, integrated over the whole numerical domain, 
for model series B5el0 (upper panel), Blell (middle panel), and 
B2ell (lower panel). Those for the MRI runs are plotted until 
the shock surface reaches the outer boundary in each model. Only 
those for model B2ellA12.5 are plotted even after the shock passes 
the outer boundary in order to see a clear saturation of the mag¬ 
netic energy of the poloidal component. The red crosses represent 
the moment when the shock surface reaches the outer boundary in 
model B2ellA12.5. 


BlellA12.5 and B2ellA12.5 are 91% and 96% of these 
values, respectively, i.e., they are close to convergence. 
Indeed, it is expected that if we were able to afford twice 
higher resolution, we could achieve convergence. 

Our results also suggest that larger initial magnetic 
fields may result in larger saturate d values. This is con - 
sistent with the results obtained by Hawley et al. (1995), 
who performed local box simulations of IVIKI in the con- 


Figure 5. Color maps for the strength of the poloidal magnetic 
fields for model B5elOA12.5 at tpb = 12 ms (upper panel) and 
200 ms (lower panel). The black line in the upper panel represents 
the shock surface. 



Figure 6. The time-averaged saturation values of the magnetic 
energy of the poloidal component (crosses) and the fitted curves 
(solid lines) with respect to the resolution for model series Blell 
and B2ell. The dotted lines represent the saturation values, ai. 


text of accretion disks. Masada et al. (2015) also claimed 
that the saturation depends on the initial fields, however, 
no resolution study was done. As shown here, the reso¬ 
lution dependence should properly taken into account in 
discussing the saturation. 

The magnetic energy of the toroidal component, 
also shows the exponential growth in each model (see Fig¬ 
ure |^. Since the non-axisymmetric MRI for the toroidal 
components cannot be treated with the current simula¬ 
tions, this is not due to the MRI but due to the winding 
of the MRI-amplified poloidal component by differential 
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runs. 

rotation. This is understood from the fact that Eb^ con¬ 
tinuously increases even after the poloidal component is 
saturated and become one order of magnitude greater 
than Eb^‘ At the end of the simulations, Eb^ still con¬ 
tinues to increase gradually in most of the models. Only 
for models BlellA25 and B fell A 12.5, Eb^ has nearly 
reached the saturated values. Incidentally, the numeri¬ 
cal convergence is achieved in Eb^ for model series Blell 
and B2ell. 


3.2. Impacts on Global Dynamics 
3.2.1. Background Runs 

With our choice of Lj^ = 1.0 erg s“^, the BG run with 
no magnetic field and rotation fails to explode, the shock 
wave being stalled at r < 150 km (see black line in Fig¬ 
ure]^. Although the shock surface is deformed by SASI- 
like oscillations during the early postbounce phase, which 
is imprinted in the zigzag evolution of the shock radius in 
Figure!^ at tpb < 300 ms, it becomes almost spherically 
symmetric later on. 

Initial rotation of |T/1F| = 0.25 % substantially 
changes the behavior of the shock evolution. In fact, flu¬ 
ids at middle to low latitudes tend to expand toward a 
larger radius thanks to centrifugal forces. The maximum 
shock radius gradually increases and exceeds 200 km by 
^pb = 700 ms (see the cyan line in Figure]^, at which 
time some parts of fluid elements are still going outward 
albeit slowly. These features are consistent with the for- 


mer findings (ISuwa et al. 2010 Nakamura et al. 

2014; 

Iwakami et ah] 2014J[ 

, that the rotation helps the explo- 

sion (bee, however, 5 

larek & Jankaj ( 

2009|). 



In the BG runs with both magnetic field and rotation, 
the shock surface propagates outward more easily com¬ 
pared with the rotation-only model, with faster propaga¬ 
tion speeds for stronger initial magnetic fields (Figure [^. 
The panels (a), (b), and (c) of Figure]^ show that the re¬ 
gions of low plasma beta (/d, the ratio of matter pressure 
to magnetic pressure) appear around the mid-latitude, 
indicating that the magnetic pressure plays an impor¬ 
tant role to push the shock outward. 

3.2.2. Dynamical Behavior of MRI Runs 

The dynamics change even more drastically when the 
spatial resolution is increased (MRI runs). Figure [^dis¬ 
plays the distributions of the plasma beta for all tne 15 
models at tpb = 578, 402, and 172 ms for model series 
B5el0 (left column), Blell (middle column), and B2ell 
(right column), respectively. In model series B2ell, all 


the MRI runs result in the formation of a collimated low- 
P jet emerging from inside the roughly-spherical shock, 
whereas the BG run yields an almost spherical expan¬ 
sion of the shock wave. Note that the low-yd region seen 
in panel (1) for model B2ellA25 at 172 ms evolves into 
a collimated jet later on as shown in panel (a) of Fig¬ 
ure Meanwhile, in a weaker-field model series Blell, 
the situation is not as simple as in model series B2ell. 
As the resolution gets higher from model Blellbg to 
BlellAlOO, the shape of shock surface changes from 
spherical to prolate, but it returns to spherical shape 
when the resolution is doubled again (see panels (b), (e), 
(h) of Figure 1^. Another doubling of the resolution, in 
turn, brings Sout the formation of a collimated low- 
yd jet emerging from inside the spherical shock (model 
BlellA25, panel (k)). In the highest resolution model 
BlellA12.5, a low-yd region is observed around the radius 
of 200 km in the vicinity of the pole (panel (n) of Fig¬ 
ure j^, which may hint at a later jet formation. Although 
the nead of low-yd region is still lingering around the ra¬ 
dius of 300 km at the end of the simulation (tpb =440 ms, 
see panel (b) of Figure]^, we expect that it would prop¬ 
agate further and eventually forms a collimated jet as 
found in model BlellA25 (see below). Finally in the 
weakest-field model series B5el0, the shock surfaces are 
roughly spherical for all the resolutions except model 
B5el0A100, which has a prolate shock, and no model 
shows a jet formation until the end of the simulation. 

We first discuss the factors responsible for the different 
shock morphorogies for different resolutions by compar¬ 
ing models BlellAlOO and BlellA25 at tpb = 213 ms, 
several milliseconds prior to the launch of the collimated- 
jet in model BlellA25. The upper panels of Fig¬ 
ure depict the distribution of the ram pressure for 
the two models. It is observed in both the models that 
a vicinity of the pole is dominated by intense downflow 
(blue region), outside of which modest outflow driven 
by relatively-low plasma beta is seen (red region; see 
also lower panels). The width of the downflow channel 
is found to be narrower for the lower resolution model 
BlellAlOO, which would be due to less effective MRI: 
the better the MRI resolved, the more efficiently angu¬ 
lar momentum is transferred outwards from the rotation 
axis, with which the rotational support decreases further 
around the pole and a broader downflow channel forms. 
Accordingly, the lower-6> edge of the outflow gets closer 
to the pole in this model, viz. the matter is ejected more 
preferentially along the pole. It is likely that this causes 
the prolate shock surface found in a model BlellAlOO. 
Note that although relatively-low plasma beta, yd ^ 1, 
is seen around the bottom of the downflow channel in 
both the models, it seems not enough to drive the mat¬ 
ter outward against the downflow (see the lower panels 
of Figure 10). The magnetically-driven mass ejection is 
only possible for the region outside the channel with such 
relatively-low plasma beta. It is found from Figure 
that the trend of broader downflow channel and thus a 
larger deflection of the outflow direction from the pole 
for higher resolution models is valid for a wide range of 
resolution in model-series B5el0 and Blell as long as a 
jet is absent. This suggest that our interpretation for the 
shock morphology is reasonable. Note that the trend dis¬ 
cussed here becomes no longer valid once a jet appears. 
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Figure 8. Color maps of the plasma beta with velocity vectors presented by arrows at tpb =578, 402, and 172 ms for model series B5el0, 
Blell, B2ell, respectively. 
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Figure 9. Color maps of the plasma beta for model B2ellA25 at tpb =195 ms (a) and model BlellA12.5 at tpb =440 ms (b), which are 
supplemental plots for panels (1) and (n) of Figure]^ respectively. 
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Figure 10. Color maps of the ram pressure (upper panels) and the plasma beta (lower panels) for model BlellAlOO (left panels) and 
BlellA25 (right panels). The ram pressure is multiplied by —1 where the radial velocity is negative. 
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since it changes the flow structure. 

According to the above discussion, the collimated-jets 
seen in some models must have been launched against 
the downflow. In the case of model BlellA25, the low- 
(3 clump around r = 70 km in the vicinity of the pole 
observed in panel (d) of Figure 10 is a prototype of the 
low-;d jet seen at a later phase y^nel (k) of Figure]^. 
We found indeed that this low-/3 clump suffers from suc¬ 
cessive depression by the downflow until it finally forms 
the collimated jet. 

To see the process of jet formation more in detail, we 
define the “low-/3 head” by the maximum radius of the 
region where the ratio of the matter pressure to mag¬ 
netic messure, each of which is angularly averaged over 
0 < is less than 0.5, and plot the evolutions of them 
for model series Blell and B2ell (see the upper pan¬ 
els of Figure 11 ). For model BlellA25 the sequence of 
the depression described above is clearly seen as an os¬ 
cillatory evolution of the low-/3 head until tpb ~ 370 ms, 
from which it grows monotonically. Model BlellA12.5 
shows a slower growth of the low-/3 head and more dis¬ 
tinct feature of oscillation, indicating that the depression 
by downflow is more significant. A similar oscillation is 
also found for model BlellAlOO, but the \ow-/3 head al¬ 
most stagnates at a small radius in this case. On the 
contrary to the weaker field case, model series B2ell 
shows no oscillation of the low-;d head, which grows al¬ 
most monotonicall y in a ll of the three models plotted in 
panel (b) of Figure |ll|f[ 

From the above discussion, the downflow seems the key 
to the propagation of low-/d head and the eventual for¬ 
mation of a collimated jet. Then the condition for the 
jet formation may be obtained by comparing the kinetic 
energy of the downflow and the magnetic energy respon¬ 
sible for the jet driving. As argued below, we indeed 
found that this comparison reasonably explains the jet 
formation. 

In the lower panels of Figure pT] we plot the time evo¬ 
lutions of the above two energie^or model series Blell 
and B2ell, defining the former as = Jy 

and the latter as = // 3 <o where the 

integrants are nonzero only for Vy < ^ and p < 0.5, 
respectively, and the integration ranges are confined to 
50 km< r < Vsh and 0 < 20°. In model series B2ell, 
the evolutions of F^Sjet appears rather similar among the 
different resolutions. In each model, F^Sjet exceeds F^df 
around tpb ~ 100 ms, which is found to approximately 
coincide with the start of the low-;d head propagation 
(see right column of Figure 11). This suggests that the 
^df-^Sjet comparison is inde^ a rough indicator for the 
jet formation. 

Unlike in model series B2ell, the values of F^Bjet 
model series Blell rather diverge among different res¬ 
olutions during a late phase after their growth nearly 
saturates around tpb ~ 250 ms. Interestingly, while re¬ 
ducing the grid size, A 50 , from 50 m to 25 m result 


The average of variable A is taken as J AdVj f dV. 

® As shown, the low-^ heads evolve more or less similarly among 
all the MRI runs of model series B2ell. Although, in the right 
column of Figure]^ all the MRI runs except for B2ellA25 show 
a trend that the jet-head radius is larger for a higher resolution 
at 172 ms, wit h w hich one may think that model B2ellA25 is an 
outlier. Figure pT] represents that there is in fact no such a trend 
on the whole time. 


in averagely-larger F^Bjet during the late phase, which 
may be simply due to smaller numerical diffusivity, an¬ 
other doubling of resolution decreases that value (see 
panel (c) of Figure |TT]). We infer that the decrease of 
^Bjet observed here is caused by the interaction of \ow-/3 
matter with downflow. Since the downflow matter in¬ 
volves high-/d (see Figure 10), the plasma beta of \ow-/3 
matter increases as it hit oy and mixed with the down¬ 
flow matter, which results in downturn of F^Bjet- Indeed, 
the panel (c) of Figure 11 shows that the downflow en¬ 
ergy, F^df 5 is averagely larger in model BlellA25 than in 
BlellA12.5, which is due to more effective MRI as dis¬ 
cussed before, and the effect of downflow is expected to 
be more standout for the latter model. Although the 
effect of the downflow basically becomes more potent 
by increasing the resolution, whether it is essential for 
the change of F^Bjet would depend on the competition 
with other factors. For the increase of F^Bjet from model 
BlellA50 to BlellA25, it is likely that the reduction of 
the numerical diffusivity is more important than the in¬ 
crement of the downflow effect. Meanwhile, the fact that 
the F’Bjet is roughly unchanged by increasing the resolu¬ 
tion in model series B2ell indicates that the downflow 
effect is insignificant in these models. One reason for 
this would be that the plasma beta of low-^d matter is 
low enough to maintain f3 < 0.5, the criterion for adding 
up F^Bjet? even after the mixing with downflow matter. 
We compare in Figure the /^-distribution of magnetic 
energy contained within^ < 20 ° for models BlellA12.5 
and B2ellA12.5 at the moment when F^df first reaches 
2.0 X 10^^ erg in each model, and found that the lat¬ 
ter model indeed involves more low-/3 matter. Another 
reason that we consider important is that models B2ell 
take shorter time to reach the saturation of magnetic 
energy than models Blell do (see Figure and lower 
panels of Figure 11 ). Since F^df gradually increases until 
attenuated by the jet formation, an early growth of mag¬ 
netic energy is advantageous to alleviate the downflow 
effect. In fact, the value of F^df the moment when it 
is caught up with by F^Bjet? generally smaller in model 


series B2ell (F^df = O-b ~ 1-2 x 10^^ erg) than in Blell 
(F;df = 1.5-2.3 X 10^8 erg). 

Bearing in mind the variations of F'Bjet ^df 

among the different resolutions mentioned for model se¬ 
ries Blell in the above, the non-monotonic dependence 
of jet formation on the resolution found in these models 
(the middle column of Figure]^ is also explained reason¬ 
ably in terms of the F^df-^Bjet comparison. For model 
BlellA50, the fact that F^Bjet almost always smaller 
than F^df is consistent with the stagnation of the low-yd 
head at small radii. Similar to model series B 2 ell, model 
BlellA25 shows the outward propagation of the low- 
P head after F^b et becomes comparable to F^df around 
tpb ~ 220 ms. Contrary to the former cases, however, 
^Bjet docs not exceed F^df so much and sometimes even 
falls behind that, as expected from the oscillatory evo¬ 
lution of the low-yd head. In the higher resolution model 
BlellA12.5, the low-yd head starts to propagate after 
F^Bjet grows comparable to F^df around tpb ~ 280 ms, 
but shows remarkable oscillations as F^Bjet occasionally 
becomes smaller than F^df by up to factor ^ 10 , due to 
the downflow effect. About 140 ms later, however, low-yd 
filaments outside the downflow choke the channel region 
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Figure 11. Evolution of low-d head (upper panels), and the downflow energy, and magnetic energy responsible for jet driving, 

(lower panels). The left and right columns are for models series Blell and B2ell, respectively. Note that the time average over the interval 
of 10 ms is taken for each plot. 
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Figure 12. /3-distribution of magnetic energy contained within 
0 < 20° for models BlellA12.5 and B2ellA12.5 at the moment 
when Edf first reaches 2.0 x 10"^® erg in each model, tpb = 198 ms 
for the former and tpb = 140 ms for the latter. The vertical-dotted 
line represents /3 = 0.5. 


(see panel (b) of Figure]^, decreasing F^df drastically and 
resulting in the acceleration of the \ow-/3 head. As men¬ 
tioned before the position of the \ow-P head is still at 
r < 300 km and no clear jet formation is observed by the 
end of the simulation. Nevertheless, since et exceeds 
F^df by almost factor 10 at that time, and the latter does 
not increase significantly afterward, we expect that a col¬ 
limated jet will form later also in model BlellA12.5. It 
should be noted that since doubling the resolution from 
model BlellA25 to BlellA12.5 renders the downflow 
effect more significant, which is disadvantageous for a jet 
formation, higher resolution runs are necessary to under¬ 
stand how the dynamics converge in terms of resolution 
for model series Blell. 

Since the jet formations discussed above take place 
close to the pole, where the coordinates become singular, 
one may be worried that the observed features are merely 


numerical artifacts. Although some level of numerical 
noises originating from the coordinate singularity may 
be inevitable in spite of the special treatment described 
in Section we believe that they are of physical origin. 
This is because the jet is born at some distance from the 
pole, ~ 10 km (panel (d) of Figure 10), which is much 
larger than the width of the region oTthe special treat¬ 
ment, and because the evolution of the low-^ region is 
reasonably understood by the above arguments. 


3 . 2 . 3 . Boost of Explosion via MRI 

Although the variation of dynamical behavior with the 
resolution seems rather complicated as described above, 
there is actually one clear trend, i.e., the faster shock 
expansion at the equator for the higher-resolutions (see 
Figure and [T^. Since the equatorial region contains a 
larger amountm mass compared to the polar region, the 
larger explosion energy is expected for the faster shock 
expansion. As shown shortly, this is indeed the case. 

Figure shows the time evolution of the diagnostic 
explosion energy, which is defined as the sum of the ki¬ 
netic, magnetic, internal, and gravitational energies over 
the fluid elements that move outward with positive en¬ 
ergies, for model series Blell. This clearly shows that 
the diagnostic explosion energy becomes larger as the 
resolution is increased. Figure [^indicates that the mag¬ 
netic effects are not necessarily lager for higher resolu¬ 
tions (e.g., compare panel (e) and (h)), which suggests 
that the magnetic pressure is not a key factor to boost 
the explosion. Note that although the collimated jets are 
driven by magnetic pressure, they give a minor contribu¬ 
tion to the explosi on energy due to their s mall volumes. 
As pointed out by Sawai & Yamada (2014), the increase 
in the explosion energy is attributes to the more efficient 
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Figure 14. Time evolutions of the diagnostic explosion energies 
for model series Blell. After the maximum shock position exceeds 
the outer radial boundary, they are plotted by dotted lines. 


neutrino heating in higher resolution models. 

How close to revival the stalled shock is roughly mea¬ 
sured by the ratio of the advection timescale, Ta, during 
which matter traverses the gain region, to the heating 
timescale, Th, within which matter gain s enough en ergy 


to overcome gravity (Thompson 2000). Following Do- 



Figure 15. Time evolutions of Ta/rh (upper panel) and the net 
heating rate per unit mass averaged over the heating region (lower 
panel) for model series Blell. See text for the definition of Ta and 
Th- 


lence et al. (2013), we define the advection timescale as 


pRgah 

Rsh 


dr 


(31) 


where the double angle bracket implies that the solid- 
angle average over Air as well as the time average over 
the interval of 10 ms are taken. is the mean shock 
radius, whereas i^gain is defined as the innermost radius 
at which the solid-angle-averaged net heating is positive. 
The heating timescale is defined as 


,nede brackets mean that 


^ -h p^)r‘^dr 

I 1 


(32) 


where the single angle brackets mean that the only solid- 
angle average is taken. The upper panel of Figure 
plots the evolution of Ta/rh for model series Blell. The 
comparison of this figure with Figure indicates that 
shock revival, which is indicated by positive explosion 
energies, roughly corresponds to Ta/rh ^ 1. It is also 
evident that higher resolutions result in higher heating 
effici ency. 


In Sawai & Yamada (2014), we argued that this is due 
to the increase of Ta in the higher resolution models as 
a result of more efficient angular momentum transfer, 
which leads to the expansion of the heating region. This 
is true of the current models. Comparison between mod¬ 
els BlellAlOO and BlellA12.5 at tpb = 180 ms shows 
that the heating region is thicker (see the upper pan¬ 
els of Figure 16) and the amount of angular momentum 
contained in the heating region is larger for the latter 
model: they are 7.0x lO^^g cm^s“^ for model BlellAlOO 
and 1.9 x lO^^g cm^s“^ for model BlellA12.5 at tpb = 
180 ms. 
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Figure 16. Top panels: color maps for the net heating rate per unit mass for models BlellAlOO (panel (a)) and BlellA12.5 (panel (b)) 
at tpb = 180 ms. Panel (c): Color map for the proton fraction, for model BlellA12.5. Panel (d): Color map for the strength of the 
poloidal magnetic fields for model BlellA12.5. 

Besides the increment of Ta, we found in this paper ^ 500): the product is YpT^{r]e) ^ 30. We 


that the reduction of Th owing to a larger heating rate 
per unit mass is also contributing to the larger Ta/rh 
in the higher resolution models. As shown in the lower 
panel of Figure 15, the heating rate per unit mass dur¬ 
ing ^ 100-250 ms, the period crucial to shock revival, 
becomes larger as the resolution increases. The com¬ 
parison of the upper panels of Figure for models 
BlellAlOO and BlellA12.5 indicates that this is orig¬ 
inated in a patch of region with large heating rates 
around the equator observed in the latter model. We 
evaluated the heating and cooling rates separately and 
found that the relatively inefficient cooling in the patch 
compared with the surroundings is responsible for the 
larger net heating rate. From Equation (|^, the cool¬ 
ing rate per unit mass, is proportional to {kT)^ 

and YnT^{—rie) + YpT^{r]e). We found that there is no 
substantial difference in {kT)^ between the patch and 
surroundings but that YnJ^^{—r]e) + YpJ'^{r]e) is several 
times smaller in the patch. In the surroundings, where 
Yn ^ Yp 0.2, J’5(-Pe) ^ 10, and 7*5(Pe) ^ 900, the 
products are YnJ^s{—r]e) ^ 2 and YpJ^^{r]e) ^ 200, viz., 
the cooling is dominated by electron capture since elec¬ 
trons are much more abundant than positrons. On the 
other hand, the electron capture is found to be relatively 
inactive in the patch due to small number of protons 
(Yp r\j 0.05, see panel (c) of Figure [l^ and electrons 


found that the positron capture rate is small as well 
in the patch: YnJ^s{—r]e) ^ 20, where Yn ^ 0.8 and 
r\j 30. To summarize, the the larger heating 
rate in the patch is caused by poverty of protons and 
electrons, and a low electron capture rate as a conse¬ 
quence. 

The low-1^ (equivalently low-Yg) region coincides with 
the location where the poloidal magnetic field is rela¬ 
tively strong (compare panels (c) and (d) of Figure 16) 


This suggests that low-1^ fluids originally locateoat 
small radii are drifted along the magnetic flux loops by 
MRI. We also found that the outflow along the rotation 
axis observed in model BlellAlOO and the collimated 
jets found in models BlellA25 and BlellA12.5 (see 
panels (e), (k), and (n) of Figure!^ respectively) also 
convey low-1^ matter from deep inside the core, which is 
reflected to the rise of the volume-averaged net heating 
rate seen after tpb ^ 350 ms for these models (see green, 
magenta, and red lines in the bottom panel of Figure 
In order to estimate the possible influences of these ef¬ 
fects on the global dynamics, we carried out two groups 
of additional test simulations based on models BlellA50 
and B2ellA50. The first one of them are extends the ra¬ 
dial outer boundaries to r = 1000 km (models BlellH 
and B2ellH). The other one is different in that the 
net neutrino heating is switched off, i.e., + Q\ 
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Figure 17. Color maps for the plasma beta for models B2ellH 
(upper panel) and B2ellNH (lower panel) at tpb = 229 ms. 

max O] (models BlellNH and B2ellNH). 

Fignre [TT] shows the profile of the plasma beta for mod¬ 
els B2ellH and B2ellNH at tpb = 229 ms, when the 
shock surface reaches r = 1000 km in model B2ellH. 
Comparing the two panels, we can immediately see the 
importance of the neutrino heating: when the shock sur¬ 
face reaches r = 1000 km in model B2ellH, that in model 
B2ellNH has just passed the radius of 500 km, and it 
arrives at r = 1000 km 65 ms later than model B2ellH. 
Even though the collimated jets observed in these two 
models are magnetically dominated, the neutrino heat¬ 
ing plays a significant role. The diagnostic explosion en¬ 
ergies at the time when the shock fronts reach the ra¬ 
dius of 1000 km are 4.6 x 10^^ and 1.8 x 10^^ erg, re¬ 
spectively, for models B2ellH and B2ellNH, implying 
that the contribution of the neutrino heating to the ex¬ 
plosion energy is about 60%. The neutrino heating is 
even more crucial for weaker initial fields. The shock 
surface in model BlellNH stays within 400 km even at 
678 ms after bounce while that of BlellH has passed 
the radius of 1000 km at tpb = 595 ms. The diagnos¬ 
tic explosion energy in model BlellH is 4.7 x 10^^ erg 
at the time when the shock front reach the radius of 
1000 km, while that in model BlellNH is negligibly 
small, ^ 2 X 10^^ erg, through the simulation. Note 
that longer time simulations following the propagation 
of the shock front through the whole progenitor would 
be necessary to correctly measure the explosion energy. 

The small contribution of magnetic field to the explo¬ 
sion energy discussed above is a substantial difference 



.y 




Figure 18. Upper panel: time evolutions of the magnetic energies 
contained in the range of 50 < (r/km) < 500 for models RinSO and 
B2ellA50. Lower panel: time evolutions of the diagnostic explo¬ 
sion energies for models RinSO and B2ellA50. After the maximum 
shock position exceeds the outer radial boundary, they are plotted 
by dotted lines. 


from previous MHD simulations involving magnetar- 
class initial fields, in which a magnetic field alone boosts 
the explosion energy up to ^ 10^^ erg, accon ipanying 


a jet with a r ather-large opening angle (e.g., Yamada 


& Sawai 2004). This implies that in our models, the 
Maxwell stress is weaker, and thus the extraction of the 
rotational energy is less efficient than in those simula¬ 
tions. 


3.3. Effects of Inner Boundary 

We have seen that the MRI efficiently amplifies weak 
seed magnetic fields to dynamically important strengths, 
having a positive impacts on explosion. In this section, 
we investigate whether the above results depends on the 
location of the inner boundary, shifting it to smaller radii 
in model B2ellA50. 

We ran a new simulation with the inner boundary at 
r = 30 km to take the effect of strong differential rotation 
beneath the radius of 50 km into account, which we refer 
to as model RinSO. The spatial resolution of this model 
is similar to that of model B2ellA50 outside the radius 
of 50 km and is 30 m at the inner boundary. The fraction 
of the volume where Ymri is less than 10 is only a few 
percent inside the radius of 50 km, which is similar to 
that outside (see Figure [^. 

The growth rates of MRI inside the radius of 50 km 
are found to be lager than those outside on average. 
The upper panel of Figure indicates that the expo¬ 
nential growth rate of Eb^ averaged over the range of 
50 < (r/km) < 500 is larger for model RinSO than for 
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model B2ellA50. This implies that the magnetic fields 
amplified inside the radius of 50 km are advected out¬ 
ward in model Rin30. Note that the degree of differen¬ 
tial rotation is even greater inside the radius of 30 km, 
and thus the above effect would be more pronounced if 
we carried out the simulation with a yet smaller radius 
of inner boundary. On the other hand, we found that 
the saturation level of Eb^ is unaffected by the position 
of the inner boundary, which may be reasonable if the 
saturation is determined by the strength of the numeri¬ 
cal diffusion, and hence the spatia l res olution, which are 
more or less the same (see Section 3.1). The evolution of 
Eb^ is similar to that of Eb^- 
The global dynamics does not change significantly, ei¬ 
ther, by moving the inner boundary position. The low- 
/3 collimated jet that emerges from inside the roughly- 
spherical shock is a feature common to both models 
B2ellA50 and Rin30. Whereas the evolutions of the jets 
are rather different between the two models, the shock 
propagation on the equator are very similar between the 
two. As shown in the lower panel of Figure the di¬ 
agnostic explosion energies of the two models are also 
nearly identical. From these results, we believe that the 
conclusions of the current study are not affected by our 
choice of the inner boundary position, r = 50 km. 


4. DISCUSSION AND CONCLUSION 

We have performed MHD simulations of the core¬ 
collapse of rapidly-rotating magnetized stars in two 
dimensions under axisymmetry, changing both the 
strength of magnetic field and the spatial resolution. Our 
goal is to study the behavior of the MRI in core-collapse 
supernovae and its impacts upon the global dynamics. 
As a result of computations we found the followings. 

The MRI greatly amplifies the seed magnetic fields 
even in the dynamical background of the core-collapse. 
Although the dominant mode, buoyant mode or Alfven 
mode, differs from location to location, the former plays 
a primary role in the exponential growth phase. It is 
true that the linear growth rate of the buoyant mode is 
independent of the wavelength, a certain degree of high 
spatial resolution seems necessary to correctly capture 
the exponential growth. The magnetic energies of the 
poloidal component gets nearly saturated within the sim¬ 
ulation times in all models, where the saturation level is 
higher for larger initial magnetic fields. The magnetic 
energies of the toroidal component grow continuously in 
most of the models, on the other hand, and the core be¬ 
comes toroidal-field dominant. 

The MRI also has a grate impacts on the global dy¬ 
namics. Models in which the MRI are well resolved show 
faster expansions of shock surface and obtain more pow¬ 
erful explosions. The formation of collimated jet is also 
found in models where the initial magnetic field is rela¬ 
tively strong and the MRI is well resolved. The following 
two effects are found to be the key to the boost of explo¬ 
sion: the first one is the expansion of the heating region 
due to the outward angular momentum transfer. This 
makes the advection timescale, or the time for matter to 
traverse the heating region, longer and thus enhance the 
heating; the second effect is the drift of low-1^ (equiva¬ 
lently low-Ye) matter along the MRI-distorted magnetic 
flux loops as well as their ejection by the jets, from deep 
inside the core to the heating region. The cooling due to 


electron capture is reduced in the low-1^ region, and the 
net heating rises as a result. 

The diagnostic explosion energies obtained in the cur¬ 
rent simulations are much smaller than the typical value 
of ^ 10^^ erg in reality. Note, however, that our choice of 
the neutrino luminosity, = 1.0 x 10^^ erg s“^, 

which is assumed to be constant, is quite modest. Ac¬ 


cor ding t o the core-collapse simulations by Bruenn et 


al. (2013), who employed the flux-limited diffusion with 


the ray-by-ray-plus approximation for neutrino trans¬ 
port, both and are ^ 5.0 x 10^^ erg s“^ at 
tpb ^ 100 ms and decay to ^ 1.0 x 10^^ erg s“^ over 
several hundred milliseconds. If such an evolution of lu¬ 
minosity is adopted in our simulations more energetic 
explosion would be obtained. 

Although our choice of the initial magnetic field 
strength, ^ 10^^ G, is much smaller than th ose as¬ 


sumed in the former global MRI simulations (e.g., Ober- 
gaulinger et al. |2006| [Shibata et al.|[2QQ6l |bawai et al, 


2013ap , most of progenitors of core-collap se supernovae 


may posses even weaker m agnetic fields (Wade & the 


MiMeS Collaboration 2014). Since lower saturation mag- 
netic helds are expected tor weaker initial magnetic fields 
according to our results, the impact of MRI in ’’normal” 
supernovae will be smaller than that found in this work. 
Although it is important to study much weaker magnetic 
fields, simulations will be computationally expensive and 
thus currently unfeasible: a reduction of the initial mag¬ 
netic field by half, with the spatial resolution kept at the 
current level, demands eight times higher computational 
cost for 2D time-explicit simulations. 

The dependence on the initial rotation also needs to 
be investigated, since the rotation spe ed of stars is also 


likely to di stribute over a wide range (Ramfrez-Agudelo 


et al. 


2013). We are currently undertaking such studies. 


and the results will be presented elsewhere in the future. 

Although our simulations are 2D under axisymmetry, 
supernovae occur in three dimensions in reality, and non- 
axisymmetric effects such as dynamo, three-dimensional 
turbulence, non-axisymmetric modes of various instabil¬ 
ities may be important. One should keep in mind that 
these effects possibly alter results obtained by current 
2D simula tions. For example, 3D-MHD simulations per¬ 
formed by Mosta et al. (2014) demonstrated that mag¬ 
netically driven-jets can be destroyed by the m = 1 
mode of a kink-type instability, whereas such destruc- 
tio n of the jet was not o bserved in 3D-MHD simulations 
by Mikami et al. (2008). In order to know how essen- 
tial non-axisymmetric effects are, 3D global simulations 
are m andatory. During t he reviewing process of this 
paper, Mosta et al. (2015) published the results of the 
first global-3D simulations of the MRI in proto-neutron 
stars. Under quadrant symmetry, they had simulated 
the evolution of the MRI for 10 ms, and found the for¬ 
mation of large-scale, strong toroidal fields, which hints 
at later magnetically-driven mass ejections. Such simu¬ 
lations have only just begun, and the possible 3D effects 
mentioned above should be studied in detail in the fu¬ 
ture. This requires long-term, large-domain, full 3D sim¬ 
ulations, which may be marginally feasible with exa-flops 
c omputers of the next gener ation. 


Masada et al. (2007) and Guilet et al. (2015) argued 
that the neutrino viscosity may hamper the growth of 
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MRI deep inside the core, i.e., r < 30 km for a fast 
rotation like ours. Applying the magnetic field of ^ 10^^- 
10^^ G obtained i n our BG run s for r < 30 km to the fast 
rotation model of Guilet et ah (2015), we found that the 
neutrino viscosity may be marginally important there. 
Since the inner boundary condition of our MRI runs is 
given by the data of the BG runs of low resolution, the 
artificial suppression of MRI by numerical diffusions may 
effectively mimic the damping by the neutrino viscosity. 
We hence believe that full-sphere simulations including 
the neutrino viscosity will not change our conclusions in 
this paper so much, if the viscous process is important 
at all. 
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